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Heat Transport in a Random Packing of Hard Spheres 
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Heat conduction in a random packing of hard spheres is studied by nonequihbrium molecular 
dynamics simulation. We find a hard-sphere random packing shows higher thermal conductivity than 
a crystalline packing with same packing fraction. Under the same pressure, the random structure 
causes reduction of thermal conductivity by only 10% from crystalline packing, which is consistent 
with the experimental fact that amorphous materials can have high thermal conductivity which is 
comparable to that of crystals. 

PACS numbers: Valid PACS appear here 



I. INTRODUCTION 

It is interesting problem how the material structure ef- 
fects on the transport property. Amorphous solid often 
shows lower conductivity than crystals. For example, the 
thermal conductivity of silica, which is a glass material 
used for window glass, is one-eighth of that of quartz, 
which is the crystal made from the same components. In 
other cases, amorphous solids of aluminum nitride show 
that the maximum thermal conductivity is 85% of that of 
pure crystal [1]. Such fact indicates the thermal conduc- 
tivity of amorphous solids can be comparable with that of 
crystalline solids. Similar feature is also found for electric 
transport in amorphous oxide semiconductors Q- 

The thermal conductivities of a low-density fluid and 
a crystal have been studied well rather than amorphous 
solid. In theoretical approaches, the linear response the- 
ory well describes the property of macroscopic thermal 
transport in crystal. Numerical approaches have been 
also done to study thermal conductivities of fluid and 
crystal from microscopic dynamics, employing molecu- 
lar dynamics technique [3, 4]. On the other hand, com- 
plexities in amorphous structures increase difficulties of 
studying thermal properties. As one of theoretical ap- 
proaches, mode-coupling theory has been developed to 
treat a localized property which comes from microscopic 
structures [5.]. However, the treatment of structural ef- 
fects is insufficient to fully understand the properties. 
Then, numerical approaches are particularly useful to un- 
derstand macroscopic properties from microscopic struc- 
tures. In particular, thermal properties of amorphous 
solids in nonequilibrium state should be revealed for ap- 
plications to novel materials. 

Before entering conduction problem, there is a long 
history of static amorphous structure. Since 1960s ^ 
ISl, structural features of jammed packings have been ex- 
tensively studied, in particular, those of the most dense 
packing called random close packing (RCP) have been 
revealed. Recently, heat transport in jammed particle 
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FIG. 1: (a) A typical packing structure obtained from our 
packing method. Colors of the particles correspond to the 
local order parameter q explained in the body, which shows 
that the structure is random, (b) The percolating cluster of 
closely- located particles. 



packings was studied to understand the dynamical prop- 
erty of glass [9^-Tl]. However, mechanisms remain unclear 
on contributions of such microscopic structures to ther- 
mal conduction in jammed state. 

The purpose of this article is to investigate thermal 
conductivities of jammed packings by the nonequilib- 
rium molecular dynamics simulations. To examine struc- 
tural effects on transport, we consider hard-sphere sys- 
tems. The hard spheres feel perfect exclusion volume 
effects of other particles, which requires special proce- 
dures to realize random packings. Since previously de- 
veloped method|12Ml4| takes not a small computational 
times, we develop two efficient methods to obtain ran- 
domly packed structures. Further, by imposing temper- 
ature gradient, we will show that thermal conductivities 
of random packings are comparable with those of crystals 
under same pressures. We also address that this is due 
to a characteristic transport paths which enhance high 
energy-transport in the random structures. 



This paper is organized as follows. In the next sec- 
tion, we improve the previous packing methods to obtain 
random packings efficiently. Then we analyze local struc- 
tures of the obtained packing in §III[ Dynamic properties 
of the packings under thermal gradient are shown in §IV[ 
The last section is devoted to summary and discussions. 



II. RANDOM PACKING 
A. Model and Methods 

In this section, we show two efficient methods to pro- 
duce random-packing. We consider hard-core particles. 
The interparticle potential is described by 



(r) 



oo r < a 
r > cr ' 



(1) 



where r and a denote a distance between particle 
pairs and the diameter of particles, respectively. Since 
random-packed structure is not a thermal equilibrium 
structure, some artificial procedures are needed to realize 
it. It is not a simple task since random putting of parti- 
cles causes overlapping which is prohibited by hard-core 
potential. To solve this some efficient packing methods 
have been developed J12l - [l^ . Some of the methods take 
the following procedures. We prepare random initial con- 
figuration of hard-core particles with sufficiently small 
diameter. Then, by increasing the diameter, we obtain 
dense-packed structures. In the original procedure, the 
diameters increase with a constant expansion rate as the 
simulation time increases. Therefore, the particles have 
spaces to expand more in each simulation step. To ex- 
pand more efficiently, we adopt two techniques explained 
below. In our simulation, we employ the event-driven 
molecular dynamics method p^. In this method, simula- 
tion steps proceed by collision events of particles, not by 
a constant time-step integration. 



1. Monodisperse packing method 

For the first method, we adopt a variable expansion 
rate to improve the original method [1^. In this method, 
the diameters of particles increase uniformly. We call it 
"Monodisperse packing method" . 

We consider L^x LyX L^ 3-dimensional box with peri- 
odic boundaries. In this box N spheres are put randomly. 
And initial radii are set to be unity. The initial veloci- 
ties of particles are randomly assigned from the Maxwell 
distribution. We adjust the total velocity of particles to 
be zero. The initial state is produced by adding each 
particle one by one into the simulation box. If a newly 
entering particle overlaps with already existing particles, 
the trial is rejected. Rejection probability is about 5%. 
Choosing the initial packing fraction to ^/(/)sc ^ 0-5, 
where ^sc denotes the volume fraction of the simple cu- 
bic packing (0sc = 7r/6 ^ 0.52). We want to expand the 



diameters of the particles keeping their random initial 
positions. To obtain dense packing, however, reconfig- 
uration is necessary because of the hard-core potential. 
This is performed by event-driven dynamics. In accor- 
dance with this scheme, we take the following steps: 

1. Find the minimum gap between the particles / over 
all the particle pairs. 

2. Increase the diameters of the particles uniformly 
by xl^ where we introduce a parameter x{< 1) as 
packing speed. 

3. Proceed the time until the earliest collision occurs. 

4. Translate all particles by half of the next collision 
time. 

5. Repeat the steps 1 to 4 until the packing fraction 
reaches an aimed value. 



2. Polydisperse packing method 

Although the above method enables us to produce 
dense-packing, it takes a lot of calculation cost. To 
solve this, we improve the packing process of the former 
method. In this second method, the expanding rates are 
different particle by particle, so we call this "Polydisperse 
packing method" . 

We first determine a final diameter crfinai- The initial 
configuration is obtained by the same way with the first 
method, and then we take following steps: 

1. Find the local minimum-gap U between the particle 

i {i = 1, ..., A^) and surrounding particles. 

2. Increase a diameter of particle a^ by xk. Here, if 
the new diameter exceeds afinai, it is set to afinai- 

3. Perform the same expansion process (steps 1 and 
2) over all particles. 

4. Repeat the steps 1 and 2 until the new diameters 
of all particles are calculated. 

5. Proceed the time until N collisions occur. 

6. Repeat the steps 1 to 4 until diameters of all par- 
ticles reach (Jfinai- 

We have introduced the packing speed x in the above 
two methods. Choosing x = 0.99 which is the maximum 
speed in our simulation, we obtain a random packing 
shown in Fig. [H This parameter controls the packing 
processes from quench to anneal explained below. 
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FIG. 2: The evolution of the packing fractions (bottom) and 
MSDs (top) as a function of the number of colhsions in the 
monodisperse packing method. The hues correspond to ^freeze 
and the maximum of MSD in the top and the bottom figures. 



FIG. 3: The frequency distribution of the order parameter q. 
The arrows at ^ = 0.92 and 1 correspond to the existence of 
the hexagonal close-packed and the face-centered cubic struc- 
tures, respectively. 



B. Packing process 

We can obtain a homogeneous packing at any simu- 
lation time by using the monodisperse packing method. 
If the aimed packing fraction is not decided before the 
simulation starts, the packing finally becomes the closest 
packing. On the other hand, when the particles jam, we 
spend a lot of calculation costs until the jammed pack- 
ing crystallizes. To observe the dynamical properties of 
the monodisperse packing method, we finish each simu- 
lation when the maximum of displacement in each cycle 
becomes less than 10~^. Figures [2] show the evolution of 
packing fraction 0/^sc and mean-square displacements 
(MSD) as a function of the number of collisions per parti- 
cle by changing x from 0.99 to 0.0001 in the monodisperse 
packing method, respectively. Here, MSD is defined as 
an average of square of the difference between the initial 
and the current positions. In each simulation the system 
size is set to L^ = Ly = 20, L^ = 40, and N = 1000. 

For the fastest monodisperse packing x = 0.99, the 
displacement of each particle is between twice and triple 
of the final radius (~ 1.34). This means that each parti- 
cle keeps almost its initial random position. In this case, 
the packing fraction approaches ^rcf/^sc = 1-22. This 
value denotes the packing fraction of randomly closed 
packing (RCP). Thus larger values of x correspond to 
the process of quench. The displacement continues to 
grow beyond 2cr as x < 0.001, which indicates structural 
relaxation toward crystallization. We note that the sat- 
uration of MSD at (L^ + Ly + I/^)/12 is the maximum 
in the present geometry, thus it does not mean freezing. 
When the expanding speeds are sufficiently slow, packing 
fraction grows beyond ^rcp and approaches the packing 
fraction of face-centered cubic (FCC). This indicates the 
system crystallizes. In this situation, packing structure 
becomes that of the closest packing, such as FCC and 



hexagonal closest packing (HCP), which are thermody- 
namically stable. In order to distinguish whether the 
packing is crystal or not, we observe a simple order pa- 
rameter for i-th particle described as 
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12^ 



cos 6 



Jh 



(2) 



where n^ denotes the number of particle in the range of 
the local gap k < O.lScTfinai- This parameter takes 1 for 
FCC and 0.92 for HCP. The angle 6^ is defined as the 
maximum bond-angle between j-th and other neighbor- 
ing bonds. Figure [3] shows a frequency distribution of q. 
For X < 0.001, the peaks at g = 1 and 0.92 prominently 
appear. In this case the packing process corresponds 
to anneal. On the other hand, the dominant peaks of 
the distribution appear below q = 0.85 as x increases. 
In this case, the packing process corresponds to quench 
since a local crystalline structure does not appear (see 
Fig.lIJa)). This crossover corresponds whether the MSD 
is larger or smaller than (2cr)^ when the packing fraction 
pass through the freezing point, 0freeze/0sc = 0.934. 

Similar crossover between anneal and quench is also 
observed in the polydisperse packing method. By us- 
ing this method, the calculation cost is reduced drasti- 
cally, thus we can obtain larger packings. Figures 2] also 
show the evolution of 0/0sc and MSD in the polydis- 
perse packing method, respectively. In each simulation 
the system size is set to L^ = Ly = 20, Lz = 400, and 
N = 10000, which is ten times longer than that used 
in former method on the z-direction. And we fix the fi- 
nal packing fraction (/)finai/</>sc = 1.203 since the packing 
fraction cannot reach to ^rcp/^sc due to the following 
reason. We note here that highly dense packing cannot 
be obtained by the polydisperse packing method. In our 
simulation, it is hard to expand all particles to final di- 
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FIG. 4: The evolution of packing fractions (bottom) and 
MSDs (top) in the polydisperse packing method. The hnes 
correspond to (/)freeze and the maximum of MSD in the top 
and the bottom figures. 



ameter for ^finai/^sc > 1.214 with x = 0.99. This is 
because some particles rapidly expand and freeze in the 
method, and therefore the rest of particles cannot suc- 
cessfully expand their diameters due to the lack of the 
free space. While the polydisperse method has huge ad- 
vantage in the calculation cost, the method is restricted 
up to the packing fraction (/>finai/</>sc = 1.214 for the 
fastest packing speed x = 0.99. 



III. ANALYSIS OF PACKING STRUCTURE 

The packings obtained by the monodisperse and the 
polydisperse packing method show the same properties 
analyzed below. Therefore, the figures in this section 
will show results only for the polydisperse case. 

A major signature of random packing structures is the 
absence of long-range order. To distinguish the packing 
structures, the radial distribution function (RDF) ^(r) is 
calculated. The function is shown for several values of 
X in Fig. [5l The packings for smaller x show crystalline 
peaks, those are expected in FCC and HCP structures. 
On the other hand, the packings for larger x show only 
two characteristic peaks at rja = >/3 and 2, which is a 
general feature of RCP[16|, [l3- We note that the peak 
at T I G = 1 corresponds to particles in contact. 

We also calculate angular distribution function P(S) 
shown in Fig. [6l where d denotes the bond-pair angle. 
We find following two characteristic features on g(r) and 
P{S) of quenched packings: While the peak at r j a — 2 
appears in g{r\ rectilinear arrangement (see Fig |71 (left)) 
seems not to exist from P(6> = tt) c^ 0. This fact seems to 
imply this peak comes from another structure, e.g. cage 
structures such as the honeycomb structures (see Fig. [3 
(right)). While the other characteristic peaks of crystals 
do not appear, the peak at rja = \fZ remains. These 



FIG. 5: Radial distribution functions for the various packing 
speeds. 



features have been also observed in previous study [17|, 
but those reasons have not been mentioned. To clarify 
the origins of those peaks, we consider two cases of con- 
tacting particles and their spatial contacting bond-angles 
illustrated in Fig. [8] and [9l 

As in Fig. [8] (top), we consider three particles. Possible 
positions for the third particle around the two contact- 
ing particles are restricted in the range 7r/3 < 6^3 < tt. 
If particles distribute homogeneously in that range, the 
probability of finding ^3, jz{Pz\ is proportional to sin ^3 
illustrated in Fig. [8] (top). In this three-particles case, we 
obtain a reweighted distribution ^3(^3) described as 



/3(^3) = 27rsin(e3). 



< 6/3 < TT 



(3) 



This reweighting also tells us why we see no peak at 6> = tt 
(see Fig. [8] (bottom)). 

Similarly, we consider four-particles unit as in Fig. [9] 
(top). We assume triangular arrangements of three parti- 
cles and then consider positions of fourth particle. Since 
the triangular arrangements corresponds to 6^3 = 7r/3 
in the three-particles case, the assumption is reasonable 
that such arrangement will be frequently observed in the 
packing. We obtain the reweighted distribution ^4(6^4) as 
the following expression: 



v^m = pmf^m-^ i < ^^ < T 







sme/4 



V'(l-cos(94)(l/2 + cos(94)* 



(4) 



As shown in Fig. [9] (bottom), the reweighted distribution 
^4(^4) is approximately fiat in the region 7r/3 < ^4 < 
27r/3 . This result indicates that the peak structure we 
see in P{S) is solely the result of the reweighting proba- 
bility in the region 7r/3 < ^4 < 27r/3 . 
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FIG. 6: Angular distribution functions for the various packing 
speeds. 
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FIG. 7: The rectihnear arrangement and the honeycomb 
structure. Both structures contribute to the peak at r/a = 2 
in g{r). 



We have found the characteristic local structure of ran- 
dom packing. In our analysis, the random packings ob- 
tained by two our methods show the same random struc- 
ture. Therefore, in the following section, we show results 
obtained only by the polydisperse packing method be- 
cause of its computational inexpensiveness. 



IV. HEAT TRANSPORT IN RANDOM 
PACKING 

A. Simulation settings 

We investigate thermal properties of the random pack- 
ing and compare to those of crystals. As the crystalline 
structure, we adopt FCC, which is the most stable struc- 
ture in three-dimensional hard-core system. Here, we 
note that the thermal conductivity highly depends on the 
packing fraction. Therefore, we use a little loose pack- 
ings both for FCC and random ones and compare them 
by setting the density to same. 

We impose periodic boundary conditions in x- and y- 
directions. And we set heat walls on the both sides of 
z-direction in the following way. At first, we make ho- 
mogeneous packings with periodic boundary condition. 
And then particles in the region z = to 2g j \fZ and 
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FIG. 8: Possible arrangement of three particles in contact 
(top). Arrow illustrates the ranges of the bond-angle. The 
reweighted angular distribution function of ^3 (bottom) . 



z = Lz — 2g I \fZ to Lz are regarded as a part of heat 
walls and those are fixed throughout the simulations. We 
choose (110) surface for FCC as heat walls. When a free 
particle collides with particles in each wall, the parti- 
cle bounces back with a new velocity randomly chosen 
from thermal equilibrium distribution of each tempera- 
ture. The velocity distribution with temperature Tb is 
described as 



/(Vn,^t,l,^t,2) = (^(Vn)V^(^t,l)V^(^t,: 
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where 'yn,^t,i denote a normal vector and orthonormal 
tangential- vectors on a colliding point of the wall particle. 
The Boltzmann constant /cb is chosen to be unity. By em- 
ploying similar heat bath, Murakami et aL[4] showed that 
the hard-core fluid systems produce the Fourier-type heat 
conduction. They also investigate the system-size depen- 
dence of thermal conductivity, which is consistent with 
theoretical predictions of the Kubo formula and long- 
time tail of autocorrelation function of heat flux. 




FIG. 9: Possible arrangement of four particles in contact 
(top). Arrow illustrates the ranges of the bond-angle. The 
reweighted angular distribution function of O4 (bottom) . 



B. Definition of physical quantities 

Our system attains a nonequilibrium steady state, 
where heat steadily flows from high-temperature to low- 
temperature sides; Th around z = and Tl (Th > Tl) 
around z = L^^ respectively In the steady state, it is 
reasonable to assume local equilibrium [4], and then tem- 
perature is defined as 



nz) 



U^' 



\ieB{z) 




(8) 



where d and v denote spatial dimension and velocity vec- 
tor of i-th particle, respectively The bracket denotes 
time average. And B{z) means a group of particles exist 
in the region, {z — 5/2^ z + S/2). 

As the same as observed in the low-density fluid case [4] 
^/^sc = 0.69, present random packings also exhibit lin- 
ear temperature profiles (Fig.fTO]). The well-scaled linear 
profile satisfies the necessary condition to evaluate the 
thermal conductivity. Heat conduction generally comes 
from energy and mass transports. In the dense packing 
state, energy transport by collisions mainly contributes 
to the thermal conductivity. Since all particles almost 
cannot move from the initial positions, the contribution 
from mass transports can be neglected. 

In the steady state, we only measure heat flux by 
the energy received from the high-temperature heat wall 
since energy conserves in the bulk of our system. Thus 
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FIG. 10: Temperature profiles in the random packing with 
different system sizes. 



the heat fiux is defined as 

n 



(9) 



k=l 



where {AE)k denotes the received energy by k-th colli- 
sion with the wall, and tn denotes the total elapsed time 
after n collisions with the wall. Thermal conductivity is 
defined as 



^{L,) 



\/T' 



(10) 



where J^ = Q{tn)/tn- 



C. Results 

Figure [11] shows the packing fraction dependence of 
thermal conductivity both for random packing and FCC. 
In this simulation, we fix the temperature (Th = 18, Tl = 



T7, 



20 for 



6) and the box size is L^ = 240, L^ 
6000 particles. Below (j)/(j)'^c = 1.15, the random pack- 
ing rapidly crystallizes since the particles have enough 
space for reconfiguration. Typically, such crystalline nu- 
cleations occur near the high-temperature wall and then 
it grows. This phenomena have been also observed ex- 
perimentally using rigid particles [i8|. Actually, crystal 
growth technique under thermal gradient is widely used 
in the field of engineering. This technique was also ap- 
plied in colloid system [19|. In the range (j)/(j)sc > 1.15, 
thermal conductivities of the random packing and FCC 
shows diverging behavior around distinct packing frac- 
tions of RCP and the closest packing, respectively. Note 
that, in random packing simulations, the system crystal- 
lizes with a small probability even near the RCP point. 

We also investigate the system size dependence of 
the thermal conductivity at a constant packing fraction 
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FIG. 11: Thermal conductivities of the random packing (RP) 
and FCC with N = 6000 against the packing fraction. 



FIG. 12: System size dependence of the thermal conductivi- 
ties of the random packing (RP) and FCC. 



DSC = 1.203, which is slightly smaller than RCP frac- 
tion, as shown in Fig. [121 We find that thermal conduc- 
tivity of random packing is far larger than FCC, which 

— 1/2 

is discussed later. In the FCC k, is proportional to Lz 
which is consistent with linear response theory, i.e. Kubo 
formula, in which thermal conductivity is described as 



tv = lim lim — - — —77 



dt'{J({))J(t')). (11) 



And autocorrelation function is supposed to have a slow 
decay t~^/^, which is commonl y c alled long-time tails 
in hard-core particle system |20|, [21|. Thus the thermal 
conductivity for finite size system shows following size 
dependence; 



n{L,) 



log L^ in 2D 
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3D 
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On the other hand, the thermal conductivity of the 
random packing decreases as Lz increases (see Fig. [12]), 
which is opposite to the above prediction. This is ex- 
plained by the interparticle distance. In the smaller sys- 
tems, it becomes more difficult to obtain dense random 
packing around RCP. Figure [13] shows obtained pack- 
ing fractions by the monodisperse packing method ex- 
plained in section [n] For N < 3000, the obtained ran- 
dom packings show lower fraction than (j)/(j)sc ^ 1.218, 
which is universally obtained for lager systems. Since the 
effective packing fraction of RCP becomes low in these 
smaller systems, the obtained random packing are effec- 
tively closer to RCP. Therefore the interparticle distance 
becomes short. On the other hand, the interparticle dis- 
tance becomes larger as system size increases. This ef- 
fective looseness causes lower conductivity than that of 
smaller systems. 

High conductivity of random packing can be explained 
by a large number of particle collisions per unit time. 
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FIG. 13: Maximum packing fractions obtained from the 
monodisperse packing method as a function of the system 
size. 



which enhance energy transport. Roughly, collision fre- 
quency of particle pair is inversely proportional to in- 
terparticle distance. Actually, the distances in random 
packing are shorter than that of FCC. Figure [HI shows 
g{r) around r/cj = 1. In FCC g{r) shows a Gaussian 
distribution around the average of the interparticle dis- 
tance. On the other hand, a diverging peak appears at 
r/cr = 1 in the random packing. This fact implies that 
such close-particle bonds form a efficient path for ther- 
mal conduction. As a consequence of this, a single large 
cluster percolating from the one side to the other in the 
z-directions is observed by assuming that neighbor parti- 
cles within the range r/cr < 1.006 belong to same cluster 
shown in Fig. [D^b). In case of FCC, a such cluster is 
only located in the low-temperature region and does not 
bridge the space between two thermal baths. 

We also investigate pressure dependence of thermal 




FIG. 14: Radial distribution functions for the random pack- 
ing (RP) and FCC around r ja — 1 with the same packing 
fraction. 



conductivities with N = 6000. The size effect of thermal 
conductivities is negligible in this condition. In Fig. {T5\ 
thermal conductivities are plotted against pressures cal- 
culated in an equilibrium state. Compared at same pres- 
sure, the thermal conductivity of FCC is higher than 
that of the random packing in contrast to the case of 
same density condition. Thermal conductivity of ran- 
dom packing is, however, reduced only by 10% of that 
of FCC. This difference is considered to come from the 
nature of the paths of sequential collisions which is rel- 
atively straight in FCC and wondering in the random 
packing. 
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FIG. 15: Pressure dependence of the thermal conductivity. 



V. CONCLUSION 



We introduce two efficient methods to obtain random 
packings. In the monodisperse packing method, we can 
produce dense random packings homogeneously. In the 
polydisperse packing method, random packings can be 
obtained quickly in the density range (j)/(j)sc < 1.214. 
Both methods have a parameter x which controls the 
packing speed. When this speed is so fast particles can- 
not diffuse sufficiently before freezing, hence the packing 
structure becomes random. The obtained structure is 
recognized by analyzing the local crystallization parame- 
ter and the radial and the angular distribution functions 
whether the packing is random or the crystalline. 

Using the obtained random packing, we investigate 
heat transport properties by imposing parallel heat walls 
at the both ends. And then we compare the thermal 
conductivity with that of the crystalline packings. Com- 
pared at the same density, the conductivity of the random 
packing is higher than that of the crystal since the perco- 
lated cluster of closely-located particles exists in the sys- 
tem. On the other hand, compared at the same pressure, 
the conductivity of the random packing is about 10% 
smaller than that of the crystal. These results suggest 
that the amorphous solids can have comparable thermal 
conductivity to crystals. The mechanism of the enhance- 
ment of the thermal conductivity by the percolated clus- 
ter is also expected in other kinds of transports, such as 
electric conduction. 
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We investigate heat transport properties of random 
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